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ABSTRACT 

Strange mode instabilities in the envelopes of massive stars lead to shock waves, which 
can oscillate on a much shorter timescale than that associated with the primary in- 
stability. The phenomenon is studied by direct numerical simulation using a, with 
respect to time, implicit Lagrangian scheme, which allows for the variation by several 
orders of magnitude of the dependent variables. The timestep for the simulation of 
the system is reduced appreciably by the shock oscillations and prevents its long term 
study. A procedure based on domain decomposition is proposed to surmount the dif- 
ficulty of vastly different timescales in various regions of the stellar envelope and thus 
to enable the desired long term simulations. Criteria for domain decomposition are 
derived and the proper treatment of the resulting inner boundaries is discussed. Tests 
of the approach are presented and its viability is demonstrated by application to a 
model for the star P Cygni. In this investigation primarily the feasibility of domain 
decomposition for the problem considered is studied. We intend to use the results as 
the basis of an extension to two dimensional simulations. 

Key words: hydrodynamics - instabilities - shock waves - stars: oscillations - stars: 
variables: other - stars: individual: P Cygni. 



1 INTRODUCTION 

SufEciontly luminous objects, such as massive stars, are 
known to suffer from strange m ode instabilities with growth 
rates in the dynamical range i Kiriakidis. Fricke fc GlatzeJ 



Il993l: iGlatzel fc Kiriakidijll993l) . The boundary of the do- 
main in the Hertzsprung-Russel diagram (HRD) above 
which all stellar models are unstable - irrespective of 
their metallicity -, coincid es with the observed Humphreys- 
Davi dson (HD) limit iHumphrevs R.M. fc Davidson K1 
I1979I) . Moreover, the range of unstable models covers 
the stellar parameters for which the LBV (luminous 
blue variable) phenomenon is observed (for a review see 
iHumphrevs R.M. fc Davidson K.I ^9^). 

The high growth rates of the instabilities indicate a con- 
nection to the observed mass loss of the corresponding ob- 
jects. To verify this suggestion, simulations of their evolution 
into the non line ar regime have been performed. In fact, for 
selecte d models iGlatzel. Kiriakidis. Chernigovskii fc Frick^ 
il999l) found the velo city amplitude to exceed the escape 
velocity (see, however. iDorfi fc Gautschvl ll2000l) ). 

To identify a possible connection between non lin- 
ear pulsations and outbursts in luminous blue variables 
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iGrott. Glatzel fc Chernigovskii (|2002|) have studied the evo- 
lution of an initial model located in the HRD well above 
the HD limit. In this study, the shocks formed in the non 
linear regime are captured by the H-ionization zone after a 
few pulsation periods. These captured shocks start to oscil- 
late rapidly with periods of the order of the sound travel 
time across the H-ionisation zone, while its mean position 
changes on the dynamical timescale of the primary, st range 
mode instability. iGrott. Glatzel fc Chernigovskii (j^TO^) have 
shown, that this shock front oscillation is of physical ori- 
gin and therefore must not be disregarded. In particular, 
the phenomenon should not be eliminated by increasing the 
artificial viscosity. We note that the representation of the 
phenomenon requires the correct treatment of extreme gra- 
dients of the dependent variables, implying their variation 
by several orders of magnitude. It is achieved by use of a, 
with respect to time implicit Lagrangian scheme. 

The rapid shock oscillations, which are confined to a 
narrow region in the vicinity of the shock front, require an 
inhibitively small timestep and thus prevent long term sim- 
ulations. In the present paper we propose an approach based 
on domain decomposition to surmount the difficulty of vastly 
different timescales in various regions of the stellar enve- 
lope and thus to enable the desired long term simulations. 
In this procedure the various domains within the envelope 
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are to be treated separately and according to their intrin- 
sic timescales. We expect this decomposition to speed up 
the calculations considerably. An even higher speedup will 
be achieved when applying domain decomposition to two 
dimensional simulations. In this sense, the present investi- 
gation may be regarded as a preliminary study for decom- 
position in two dimensions. 

The basic equations and assumptions are introduced in 
Section 121 The domain decomposition approach is discussed 
in detail in Section |3 including a derivation of criteria for 
domain decomposition and an investigation of the proper 
treatment of the resulting inner boundaries. Moreover, tests 
of the approach are presented there. Its viability is demon- 
strated by application to a model of the star P Cygni in 
section 2] Our conclusions follow. 



2 BASIC EQUATIONS AND ASSUMPTIONS 

The evolution of instabilities of a stellar envelope is followed 
into the non-linear regime assuming spherical symmetry and 
adopting a Lagrangian description, i.e. the independent vari- 
ables are the time t and the mass Mr inside a sphere of ra- 
dius r. The equations to be solved (see, e.g.. iCoa Il980l) ) are 
given by mass conservation, 

= (1) 

momentum conservation, 

9V 2 dp ^ GMr 
+ 47rr — — — H — 



dMr 



= 



energy conservation, 

dL p dp dE _ 

dM^ ^ 'i^'dt ^ 1h 

and the diffusion equation for energy transport, 

dT 3k,{L ~ Lconv) 



dMr 



= 



(2) 



(3) 



(4) 



where p, p, T, L, and E denote density, pressure, tempera- 
ture, luminosity and specific internal energy, respectively, a 
is the Stefan-Boltzmann constant, G the gravitational con- 
stant and c the speed of light. We emphasise that ^ is 
the substantial time derivative. F or the opacities k, the lat- 
est v e rsions of the OPAL tab les ilglesias. Rogers &: Wilsoiil 
I1992I : iRogers fc Iglesiaj Il992h have been used. Convection 
is treated in the standard frozen in approximation (see 
iBaker fc KippenhahnI (|^6^), i.e. the convective luminos- 
tity Lconv is kept constant and equal to the luminosity of 
the hydrostatic initial model. Since the instabilities are lo- 
calized in the outer envelope, the evolution of the core can be 
neglected. Its properties are taken into account by imposing 
time independent boundary conditions (e.g. by prescribing 
luminosity, [vanishing] velocity 11 and [constant] radius) at 
the bottom of the envelope. For the envelope model, the nu- 
clear energy generation rate e vanishes. At the outer bound- 
ary, the gradient of heat sources is required to vanish (_F is 
the heat flux): 



grad(div_F) — 



(5) 



It is chosen to ensure that outgoing shocks pass through 
the boundary without reflection. The set of boundary con- 
ditions prescribing values for the velocity and luminosity at 
the bottom and values for the pressure and temperature at 
the top of the envelope will be denoted by {v, L) {p, T) in the 
following. 

The numerical code relies on a Lagrangian, with re- 
spect to tim e implic i t, full y co nservative difference scheme 
proposed bv lFralevI il968l) and lSamarskii fc Popovl il969l) . 
As a consistent extension the two dimensional version of 
the code currently under development (also implicit and 
Lagrangian) is based on th e method of support o pera- 
tors originally sugges t ed by ISamarskii et. alj (|^^ and 
lArdelvan fc GushchinI il982l) . For a thorough d escription 
of the method of support operators we refer to IShashkovl 
il996l) . To handle shock waves, artificial viscosity is used. 
Concern ing tests of the code, we adopted the same cr i- 
teria as iGlatzel. Kiriakidis. Chernigovskii fc Frickel il999l) . 
We emphasise, that conservativity of the numerical scheme 
is of fundamental importance when simulating instabilities 
in a stellar envelope. Considering the distribution of inter- 
nal, gravitational and kinetic energy we find, that the kinetic 
energy can be smaller than, e.g., the gravitational energy 
by several orders of magnitude. Appropriate simulations of 
stellar instabilities require the correct representation of the 
kinetic energy, and therefore energy conservation with high 
accuracy is indispensable. The difference scheme adopted 
here guarantees energy conservation. 



3 DOMAIN DECOMPOSITION 

3.1 Motivation for domain decomposition 

The stellar envelope model representing a massive star with 
mass M = SOM©, effective temperature T^s = IQ'^K and 
luminosity L = 1.17 ■ 10^L(7 ) , whic h was considered by 



This boundary condition implies (by using eauations lll3l and 
|1J| boundary values for the temperature T and pressure p. 



iGrott. Glatzel fc Chernigovskii (|2003), suffers from strange- 
mode-instabilities. These cause pulsations with velocity am- 
plitudes of 0.5 Dcsc and inflate the envelope to 2.5 initial 
radii. After several pulsation periods a shock front is cap- 
tured in the H-ionization zone. The latter is prone to sec- 
ondary instabilities and oscillates on very short timescales 
connected to the sound travel time across the front. These 
instabilities a re caused by the stratiflcation, but n ot driven 
by buoyancy iGrott. Glatzel fc Chernigovskil2002ll . Resolv- 
ing the rapid oscillations of the shock front reduces the 
timestep to very small values. This is computationally ex- 
tremely expensive and effectively inhibits the desired long 
term study of the system. 

In order to enable the treatment of the problem the 
integration interval is decomposed such that the small and 
quickly varying shock region is integrated with small time 
steps and the remaining major part of the envelope is calcu- 
lated using large time steps. This strategy of domain decom- 
position is com mo n in computationa l fluid dynamics (see, 
e.g.,|^l(l99i) and lWu fc ZoultoOOh ). For the treatment of 
interfaces between the various domains and the associated 
inner boundary conditions for purely hyperbolic systems, 
we refer the reader to these publications. We stress, how- 
ever, some fundamental differences to the previous studies: 
One of them concerns the different character of the system 



Simulation of stellar instabilities using domain decomposition 3 



1 1 


1 




1 1 


■ ■ — 


^1 







100 



-200 



D 



1400 



600 
Giridpoimis 




400 oOO 
Giridpoimis 



Figure 1. Normalized time derivatives of the velocity v as a function of gridpoint, for two typical states of the system. Decomposition 
into three (a) and two (b) subdomains (D;, X'sh, fr) as indicated by the dashed lines is suggested. 



of equations. So far, only purely hyperbolic systems have 
been considered, whereas we apply domain decomposition to 
a composite system of hyperbolic and parabolic equations. 
Moreover, we adopt a numerical scheme which is implicit in 
time. Again, this is in contrast to previous studies which use 
explicit schemes. 



the shock region. Accordingly, the zones below and above 
the shock are determined by 



T>i = {n : n < min{Vsh)} 
Dr — {n : n ^ max(I>sh)} . 



(7) 
(8) 



3.2 Criterion for domain decomposition 

In this section a criterion for the proper choice of the bound- 
aries of the various domains evolving on different timescales 
and therefore treated with different time steps will be pre- 
sented. Considering the time derivatives of various physical 
quantities the velocity i; is found to vary most rapidly and 
therefore determines the time step. 

Figure shows the time derivatives of the velocity for 
two typical states of the system, where Figure 0a corre- 
sponds to rapid shock oscillations around gridpoint 310. In 
this case, domain decomposition as indicated by the dashed 
hues is suggested. Figure 0b represents a situation, in which 
the outer envelope is collapsing onto the shock (at gridpoint 
410). Rapid variations are now also found above the shock. 
Accordingly, decomposition into two domains (as indicated 
by the dashed line) seems appropriate. Depending on the 
state of the system, we therefore need to split the domain 
of integration into two or three subdomains. 

The size of the various domains is determined by com- 
paring the time derivatives of velocity and temperature with 
the corresponding derivatives on the shock front. Therefore, 
as a first step, the position of the shock front has to be de- 
termined. For the models studied, the latter is defined by 
the maximum temperature gradient. The boundaries of the 
shock zone are then defined by the requirement, that the 
time derivative at their position corresponds to a given frac- 
tion 1/fc of the time derivative at the shock. In other words, 
the set of gridpoints Vsh belonging to the shock zone may 
be characterized by 



Vsh = 



1 



dv 
dt 



A 



dT 



dt 



> 



k 



dT 



dt 



(6) 



where n denotes the number of the gridpoint and derivatives 
with index sh correspond to their maximum values found in 



The parameter k corresponds to the ratio of timesteps in 
Vi^r and I>ah, which favours large values of k. On the other 
hand, the size of I?sh increases with k, suggesting smaller 
values. Therefore optimum results will be obtained for a 
mean choice. For the considered model a value of fc = 15 
turned out to be satisfactory. If the size of the domains T?i 
or Dr drops below a given value, they are considered to be 
part of the shock zone, and we arrive at a decomposition into 
two subdomains. For the model considered the whole grid 
consists of 512 gridpoints and the zones have a minimum 
size of 64 gridpoints. 

After decomposition, the quickly varying region "Dsh 
is integrated first with timesteps ri , r2 , . . . , r„ . Then the 
domains Vi and /or Vr are integrated with the timestep 
T — Ti + T2 + ■ ■ ■ + T„. The decomposition implies artificial 
boundaries and boundary conditions between the domains 
Vi and T>sh and T>r and Vsh, respectively. An inconsistency 
is introduced, if in a first approach the explicit inner bound- 
ary conditions are either kept fixed, or linearly extrapolated 
in time. Both cases are contained in the following extrapo- 
lation presciption: 



now = '^°"^ r' + y ) +(i-ay)-y (9) 

Told 



Y stands for p, T, v, L, and ay £ [0, 1] is a free extrap- 
olation parameter, which may be chosen independently for 
each variable. Y refers to the current, Yoid to the previous 
value of the variable. Tom is the last, r' the current timestep. 
How the inconsistency may be treated will be discussed in 
the following sections. 

Once the integration of all domains has been performed, 
one subsequent timestep is done without decomposition. We 
shall refer to it as the relaxation timestep. Relaxation pre- 
vents accumulation of residual errors. 
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3.3 The iteration procedure 

A mathematically consistent way of integrating the different 
domains, which solves the problem of artificial fixed bound- 
ary conditions is the following iterative procedure: 

(i) The computation is started at time t and Di , Dsh and 
Dr are integrated with fixed boundary values {v, L)t and 
{p,T)t for each domain. The subscript t denotes values at 
time t. 

(ii) New boundary values {v, L)t+T and (p, T)t+T are ob- 
tained by integration of the adjacent domain. With these 
boundary values the integration |(i)| is repeated. 

This iterative procedure implies implicit boundary condi- 
tions and sucessively eliminates the errors introduced by 
the artificial boundaries. 4-5 iteration cycles produce sat- 
isfactory results. However, this approach is computationally 
even more expensive than integrating the entire domain with 
small timesteps ti, . . . , r„. Concerning the present problem, 
it is therefore not relevant and has been applied only for 
comparison with other methods discussed below. Ifowever, 
for the corresponding two dimensional problem the compu- 
tational effort might be significantly reduced by the pro- 
cedure as the inversion of the most ill-conditioned matrix 
occurring there is an A^'^-process. 



3.4 Overlapping domains 

One method to reduce the error caused by the fixed inner 
boundary conditions consists of using overlapping domains 
of computation. We emphasise that the system considered 
here is composed of hyperbolic and parabolic differential 
equations. We first consider the hyperbolic part, i.e., the 
mechanical equations. Any errors or perturbations in this 
part produced at the boundaries of the domains propagate 
with finite speed into the domain of computation along the 
characteristics of the equation. For our set of equations, per- 
turbations propagate with velocities «±Cs, where Cs denotes 
the speed of sound. If the domains overlap such that the time 
for error propagation across the overlap is larger than the 
integration timestep r, we may discard the flawed values 
close to the boundary and keep only the correct values for 
the subsequent timestep. This condition implies 

where the sum of individual sound travel times is to be taken 
over all overlapping cells, h denotes the thickness of cell i. 

Concerning the parabolic part of equations I1I4I i.e. the 
diffusion equation for energy transport, perturbations travel 
across the grid with infinite speed. Therefore, the procedure 
suggested for the hyperbolic part of the equations can in 
principle not be carried over to the parabolic part. Rather 
physical quantities change on the diffusion timescale, which 
is given by Tdiff — 1^ pCp-^^ {cp is the specific heat at 
constant pressure). Consequently, we expect the effect of 
the perturbations to be small far from the boundary if the 
timestep r is sufficiently small, i.e. : 

T <^hp,Cp^j-^ (11) 



where the sum of individual diffusion times again extends 
over all overlapping cells. 

In Figure 121 a snapshot of the sound travel time (Figure 
|21a) and diffusion time (Figure |21b) across a cell is given 
as a function of gridpoint. The sound travel time across the 
overlapping region is larger than a typical timestep (approx- 
imately 1 per cent of the global free fall time) for an overlap 
of ~ 8 cells (condition llUH . However, while the diffusion time 
across the overlapping region is bigger than the timestep 
even for a small overlap in the bottom part of the envelope, 
condition II II cannot be satisfied in the outer envelope for a 
reasonable number of overlapping cells. The latter is due to 
the small heat capacity there (implying the ratio of thermal 
and dynamical timescales to be small) and in accordance 
with the validity of the non-adiabatic-reversible approxima- 
tion (NAR-approximation), which has been shown for this 
partic ular stellar model by Grott, Glatzcl & ChorniEovsk| 
l|2002f) . How satisfactory results may be obtained even if 
condition II II is not satisfied will be discussed in the follow- 
ing section. 

3.5 Inner boundary conditions 

On the basis of the domain decomposition procedure de- 
scribed above, various inner boundary conditions and their 
consequences will be investigated in this section. In any case, 
overlapping domains have been used. The tests presented 
here have been performed at various times with similar re- 
sults. Therefore the results may be regarded to be indepen- 
dent of the particular state of the system. In the tests, the 
luminosity L turned out to be the most sensitive quantity 
concerning the errors introduced by the inner boundaries. 
This is due to the fact, that in our difference scheme it is 
only accurate to first order. Therefore L is required to be 
reproduced satisfactorily compared to the results of the ap- 
proach without decomposition. 

In Figure|3the luminosity is given as a function of grid- 
point around the boundary between T^i and ©sh, which is 
more sensitive than the boundary between Tish and Dr with 
respect to the decomposition procedure. For comparison, 
results obtained without decomposition are shown as solid 
lines. Dotted lines correspond to results with decomposition, 
where the left and right columns illustrate those before and 
after relaxation, respectively. 

(i) Figure |21a shows the result obtained using the 
{v, L){p,T) boundary condition, i.e., velocity and luminos- 
ity were prescribed at the left, and temperature and pressure 
at the right boundary, respectively. This condition implies 
a discontinuity of the luminosity (al) and leads to an unac- 
ceptable relative error (~ 15 per cent) after relaxation (a2). 

(ii) Figure l^b shows the result obtained using the 
(v, L)(v, L) boundary condition. The latter is motivated by 
the discontinuity of the luminosity for the previously dis- 
cussed {v, L){p,T) boundary condition. Rather than the 
considerable discontinuity of L we expect a more tolerable 
discontinuity of its derivative for the {v, L){v, L) boundary 
condition. In fact, the results given in figure |21bl are satis- 
factory. After relaxation the relative error in the luminosity 
is of the order of 10~^ (figure |H]b2). However, in this case 
the error of L directly enters the boundary conditions for 
the subsequent timestep and leads to an accumulation of 
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Figure 2. Snapshot of the sound travel time Tcour (a) and the diffusion time r^iff (b) across a cell as a function of gridpoint. 



the error at the interface. This problem can be removed 
partially by switching between different interfaces for subse- 
quent timesteps. The remaining error can then spread sufB- 
cently and does not influence the further integration signif- 
icantly. 

(iii) Figure |21c shows the result obtained using the itera- 
tion procedure described in section HOI After four iteration 
cycles it is comparable to that obtained with the {v, L)(v, L) 
boundary condition. By performing more iteration cycles the 
accuracy could be improved even more. However, the iter- 
ation procedure is computationally much more expensive 
than the alternatives discussed and therefore of no practical 
use. 

We thus conclude, that using the domain decompo- 
sition procedure together with overlapping domains and 
{v, L)(v, L) boundary conditions, and switching between dif- 
ferent interfaces, yields satisfactory results at low computa- 
tional cost. 



3.6 Validation of the domain decomposition 
procedure 

The domain decomposition procedure with overlapping do- 
mains and {v, L){v, L) boundary conditions has been com- 
pared for validation with the original approach using no de- 
composition. The comparison starts at 6.47 • 10^ sec, i.e., 
well after the formation of the shock and the associated in- 
stability, and extends to 6.58 ■ 10^ sec. The typical periods 
of unstable strange modes driving the pulsations of the star 
are of the order of 5 • 10^ sec, whereas the modes carrying 
the stratification instabilities of the shock front have periods 
of ~ 10^ sec. Therefore, the test covers approximately 0.2 
periods of the overall envelope pulsations and 10 periods of 
shock oscillations. 

Convergence and error control is done using the follow- 
ing criterion based on a Z^-norm: 

^(E/^)^^ (12) 

where fi denotes the relative error of a physical quantity in 
gridpoint or cell i and A'^ the total number of gridpoints. E 
is the prescribed error bound and the sum extends over all 
gridpoints. contains the weight-function of the i^-norm 



which is chosen to be proportional to the mass of the cor- 
responding cell, fi (X rrii. Thus the various regions of the 
star contribute to the error in a different way and domain 
decomposition using the same error bound in all domains 
will result in different accuracies for the various domains 
and compared with the approach without domain decompo- 
sition. Accordingly the error bound has to be adapted to the 
domains contribution to the numerical error. On the other 
hand, as an identical error control cannot be guaranteed, re- 
sults with and without decomposition are expected to differ 
slightly. 

In the following figures 2] and |S] solid lines correspond 
to the approach without decompostion, dotted lines refer to 
decomposition. Figure|l|shows the velocity u as a function of 
gridpoint after 5-10^ sec (figure|31al) and 1.1-10® sec (figure 
2]bl) of simulated time, respectively. For comparison, the ve- 
locity and the density after 1.1 ■ 10® sec of simulated time are 
also presented as a function of relative radius in figures (b3) 
and (c), respectively (without decomposition). The shock is 
located at r/R « 0.95 and resolved by approximately 150 
gridpoints. The high resoluti on of the shock zone is necessary 
to rep resent its oscillations llGrott. Glatzel fc Chernigovskil 
J2002M . 

Considering 5 • 10^ sec of simulated time, domain de- 
composition yields excellent agreement with the original 
approach up to the position of the shock front (at grid 
point 320). Around the shock front the results differ slightly, 
whereas above it the agreement between the two approaches 
is again satisfying. The agreement may be found to be even 
better, if a phaseshift in time of about 10* sec between the 
results discussed is taken into account, i.e., if results of the 
original approach are compared with those obtained lO** sec 
later with decompostion (figure 2]a2). We emphasise, that 
the time interval corresponds to five oscillation cycles of the 
shock front instability, and implies a phaseshift of only 
With respect to the fact that the phaseshift is physically 
irrelevant, we thus regard the agreement as fully satisfying. 

After 1.1 ■ 10® sec of simulated time the differences be- 
come more pronounced, in particular in the vicinity of the 
shock (now at gridpoint 300). Similar to the previous dis- 
cussion, however, including a suitable phaseshift of 5 • 10* 
sec a reasonable agreement may be achieved (Figure 2]b2). 

Figure El shows the gravitational (Figure El a), kinetic 
(Figure Elb), thermal (Figure Elc) and total energy (Figure 
Eld) of the envelope as a function of time. We note that the 



6 M. Grott et al. 




270 275 

Gridpoints 



77 


4.42 


(U 










4.4 












4.38 


o 






J 


4.36 




4.34 




4.32 




4.3 




4.28 




4.26 



(a2) 



255 



260 265 



270 275 

Gridpoints 




275 

Gridpoints 



o 4.42 

c« 

"So 44 
4.38 

J 4.36 
4.34 
4.32 
4.3 
4.28 
4.26 



(b2) 



255 



260 265 



270 275 

Gridpoints 




255 



260 



265 



270 



275 

Gridpoints 



255 



260 



265 



270 275 

Gridpoints 



Figure 3. The luminosity as a function of gridpoint around tfie boundary between 15; and at gridpoint 264. Results obtained 

without decomposition are shown as solid lines. Dotted lines in the left and right columns correspond to decomposition without (1) 
and with (2) relaxation, respectively, {v, L){p,T) and (v, L){v, L) boundary conditions have been used in (a) and (b), respectively, the 
iteration procedure (4 cycles) has been applied in (c) with the {v, L)(p,T) boundary condition. 



kinetic energy given in figure |S]b is seven orders of magni- 
tude smaller than either the gravitational or thermal energy. 
It may be referred to as the energy of the pulsation and is 
therefore of central interest in the present context. Its varia- 
tion with time is one order of magnitude smaller than that of 
the gravitational energy, which itself is one order of magni- 
tude smaller than that of the thermal energy. Therefore the 
variation of the total energy is dominated by the thermal 
energy. 

Domain decomposition reproduces the gravitational en- 
ergy perfectly. This can be attributed to the use of overlap- 



ping domains, since without them considerable disagreement 
is found. The latter consists of oscillations of the solution 
around the curve given in Figure |^a. With respect to ki- 
netic (Figure |S]b), thermal (|S]c) and total energy (I5]d), the 
agreement is satisfactory up to the time 6.55 • 10^ sec. The 
deviations at later times can be partially attributed to the 
phase shift discussed above. 

To summarise, our comparisons prove the domain de- 
composition procedure to provide reliable and satisfac- 
tory results. We emphasise that domain decomposition vio- 
lates the conservativity otherwise inherent in the numerical 
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Figure 4. The velocity as a function of gridpoint with (dashed Une) and without (soUd line) domain decomposition after 5 ■ 10® sec 
(al) and 1.1 • 10^' sec (bl) of simulated time. All simulations start at t = 6.47 ■ 10^ sec (after the formation of the shock front) with 
the same initial model. In Figures (a2) and (b2) velocities obtained with decomposition arc taken 10^ sec (a2) and 5 • 10"* sec (b2) later 
than their counterparts without decomposition, respectively. This physically irrelevant phaseshift reduces the differences between the 



two approaches significantly. For comparison, the velocity and the density after 1.1 ■ 
of relative radius in figures (b3) and (c), respectively (without decomposition). 



10 sec of simulated time are presented as a function 



scheme. How this violation of conservativity contributes to 
the discrepancies discussed is an open question. 



3.7 Speed-Up of the calculations 

In order to estimate the speed-up achieved by the domain 
decomposition procedure, we assume that the number of it- 
erations I needed to solve the implicit equations with pre- 
scribed accuracy is not changed by decomposition, i.e., for 
convergence I iterations are required in the shock region Vgh 
integrated with timestep t and the same number of iterar 
tions / are needed to integrate the domains Vi and Vr with 



time step k ■ t. Moreover, we assume that / iterations arc 
required to integrate with timestep r and without decom- 
position. 

We denote the number of gridpoints by N and the size 
of the domains "Di, T>sh and T>r by A'^;, Nsh and Nr, respec- 
tively. (Ni, Nah and A',, include the overlap.) Integrating k 
timesteps r without decomposition then requires 

Oi = k-N-I (13) 

operations. Decomposing the grid into two or three domains. 
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m''^^] tfiO^secl 

Figure 5. Gravitational (a), kinetic (b), thermal (c) and total (d) energy as a function of time. Solid and dashed lines correspond to 
simulations without and with domain decomposition, respectively. 



we need for the integration of the same time interval 

02 = Ni- I + k- N,h- 1 and (14) 

03 = {Ni + Nr)-I + k-Nsh-I (15) 

operations, respectively. Thus we expect a speed-up of the 
calculations by decomposition by a factor of 



S2 



S3 = 



Oi 

O2 
Oi 

03 



kN 



-kNsh 
kN 



and 



(iV, +Nr)+ kNsh ' 



(16) 
(17) 



respectively. It essentially depends on Nsh and k. For large 
k it is approximately given by N/Nsh- However, for large k, 
Nsh increases too (see section T^. 21 . A comparison between 
the estimated and measured speed-up is presented in Ta- 
ble For the hypothetical case of large k, Nsh ~ 16 and 
an overlap of 8 gridpoints we obtain a speed-up by a fac- 
tor of 16. This situation could in principle be realised for 
shock oscillations with smooth spatial structure which can 
be represented by a small number of gridpoints. This exam- 
ple demonstrates the power of the method when applied to 
a suitable situation. 

Compared to the one dimensional case, a considerably 
higher speed-up is expected if decomposition is applied to 
two dimensional problems, since it reduces the size of ma- 
trices to be inverted, the latter being a A'^'^ operation in the 
worst case for iterative methods and the matrices consid- 
ered. 



Domains 


Size of Vsh 


Overlap 


Estimated 
Speed-Up 


Measured 
Speed-Up 


3 


160 


16 


2.75 


2.42 


3 


160 


32 


2.72 


2.36 


3 


192 


16 


2.34 


2.15 


3 


192 


32 


2.32 


2.14 


2 


272 


16 


1.77 


1.61 


2 


272 


32 


1.76 


1.56 


2 


288 


16 


1.68 


1.67 


2 


288 


32 


1.67 


1.61 



Table 1. Estimated and measured speed-up for a grid with A' = 
512 points and decomposition into 2 and 3 domains with two 
different overlaps. 



4 APPLICATION TO A P CYGNI MODEL 

We apply the method discussed above to a stellar en- 
velope model with parameters close to those observed 
for the star P Cygni. Concerning luminosity, effec- 
tive temperature and chemical composition for thi s ob- 
ject, various auth o rs llNaiarro. Hillier fc Stahll (^93) and 
iPauldrach fc Pulj Jl99CD l agree that these parameters 
should lie in the vicinity of L = 752.5 • 10^ Lq , T^a = igSOOiC 
and X = 0.31, Y = 0.67, Z = 0.02. The most un- 
certain parameter of P Cygni is its mass. Standard stel- 
lar evolution calculations in dicate a mass of M = 50Mq 
Eid fc Hartmanin (^9^), whereas spectroscopic obser- 
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Figure 6. Density p as a function of radius (in units of the initial radius) and time for an envelope model of P Cygni (top panel). 
The corresponding contour plot is given in the bottom panel. Strange mode instabilities, but no shock oscillations are resolved in these 
diagrams. 
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Figure 7. Same as Figure|n| but for the temperature T. 
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Figure 8. Luminosity in units of the initial luminosity as a function of radius (in units of the initial radius) before the shock capturing. 
The corresponding contour plot is given in the bottom panel. Note that the luminosity variations take place on the dynamical timescale. 




Figure 9. Same as Figure|H]in the vicinity of the captured shoclc and for a time interval appropriate to resolve the shock oscillations. 



Simulation of stellar instabilities using domain decomposition 13 




Figure 10. Luminosity in units of the initial luminosity as a function of radius (in units of the initial radius) and a time interval 
appropriate to resolve the shock oscillations (top panel). The corresponding contour plot is given in the bottom panel. Note that the 
rapid photospheric luminosity perturbations are defined at the captured shock (at r/Ro ~ 1.5). Below the shock, luminosity perturbations 
are governed by strange mode instabilities operating on the dynamical timescale. 
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vations are consist e nt wi th masses down to AI — 23Mq 
JPauldrach fc Pulj Jl99Ci) ). For our simulations we adopt 
M = 26.5M0, a value supported by the observation. Even 
with a more conservative mass of M — 50Mq, the model 
is known to be unstable w i th res pect to strange modes 
l|Kiriakidis. Fricke fc GlatzeJ jlOQ.^T l'l. The higher value of 
L/M adopted here amplifies the tendency towards insta- 
bility through shorter thermal timescales. 

The simulation of strange mode instabilities in P Cygni 
starts from the envelope model in hydrostatic equilibrium 
without any external perturbations. Strange mode instabil- 
ities develop from numerical noise, pass the linear regime 
of exponential growth and form multiple shocks in the non- 
linear domain. One of these shocks is captured in the H- 
ionisation zone and starts oscillating on timescales of the 
order of the sound travel time across the front (~ 0.5 days), 
whereas its mean position varies on the dynamical timescale 
(~ 10 days). In this phase of the evolution on two different 
timescales domain decomposition is used to speed up the 
calculation considerably. 

Figure |S| shows the density p as a function of radius 
(in units of the initial radius) and time for the envelope 
model of P Cygni. The corresponding contour plot is given 
in the bottom panel. Note, that contour lines here and in the 
following are not always closed, since during the evolution 
of the star e.g. the density drops to values lower than those 
of the initial model. After reaching the non-linear regime 
at t —'^ 20 days, shocks are formed in the outer envelope, 
travelling outwards and inflating the envelope successively 
up to 4.5 initial radii. During this period, the surface velocity 
reaches 55 per cent of the escape velocity at t = 75 days. 
After ~ 80 days the envelope starts to coUaps, and a shock 
originating at r/i? ~ 0.6 at time t = 70 days is then captured 
in the H-ionisation zone around r/R ~ 1.5 at f « 85 days. 
Subsequent shocks, generated by the primary strange mode 
instability, are confined to the region below the captured 
shock. 

FigureQis the analogue to Figure^lfor the temperature 
T. The shocks discussed above can easily be identified in the 
contour plot. Note the steep temperature gradient at the 
captured shock after its formation. 

Figure |H| shows the luminosity in units of the initial 
luminosity as a function of radius (in units of the initial 
radius) and a time interval before the shock capturing. The 
luminosity varies on the dynamical timescale by 10 per cent, 
corresponding to a bolometric variation of ~ 0.1™. It is de- 
fined in the inner envelope and remains constant with re- 
spect to radius above. There, luminosity perturbations can- 
not be sustained due to the low specific heat and the as- 
sociated short thermal timescales. These luminosity varia- 
tions re mind of the microvariations in P Cygni, d escribed 
by, e.g.. Ide Groot. Sterken fc van GenderenI (|200J) . These 
authors report on cyclic behaviour of the visual brightness 
with amplitudes of ~ 0.1™ and cycle lengths between 10 and 
25 days, best fitting a quasi-period of 17. S'*. 

Figure |5| shows the density p as a function of radius (in 
units of the initial radius) in the vicinity of the captured 
shock and a time appropriate to resolve the shock oscilla- 
tions. In particular, the corresponding contour plot in the 
bottom panel shows variation on two different timescales. 
It exhibits both shock oscillations with a mean period of 
~0.5 days and the evolution on the dynamical timescale. 



as indicated by the variation of the shock position between 
r/Ro « 1.56 and r/Ro ~ 1.54. 

Figure 1101 shows the luminosity in units of the initial 
luminosity as a function of radius (in units of the initial ra- 
dius) and a time interval apropriate to resolve the shock os- 
cillations. The rapid photospheric luminosity perturbations 
are defined at the captured shock (at r/Ro ~ 1.5) and re- 
main constant above it (see discussion of Figure (HJ. The 
variation amounts to 20 per cent, corresponding to ~ 0.2™ 
bolometrically. Below the shock, luminosity perturbations 
are governed by strange mode instabilities operating on the 
dynamical timescale. The photometric luminosity perturba- 
tion is therefore a superposition of two effects, dominated 
by the fast oscillations induced by the shock. We empha- 
sise, however, that the one dimensional calculations pre- 
sented here have to be interpreted with caution, since mas- 
sive sta rs are known to su f fer a lso from non-radial insta- 
bilities JClatzel fc Mehren I liggdl l. Whether the captured 
shocks survive the deformations induced by non-radial in- 
stabilities or become unstable themselves with respect to 
non-radial perturbations, remains to be tested by at least 
two dimensional calculations. They might break apart and 
could then contribute to the entire variation by stochasti- 
cally adding high-frequency perturbations to the cyclic per- 
turbations on the dynamical timescale induced by strange 
mode instabi lities. On the other hand, fr om the observations 
discussed bv lde Groot. Sterken fc van Genderen ( 20Q1) , no 
indications for the stability properties of the captured shock 
found here can be inferred, since the time resolution of the 
data is not sufficiently high (~ one measurement per day). 

Another effect of the primary strange mode instability 
consists of a mass transfer from the instability region into the 
outer parts of the stellar envelope. Owing to the Lagrangian 
approach chosen here, this implies a reduced spatial res- 
olution in the inner part of the stellar envelope. Simulta- 
neously, grid points are concentrated around the captured 
shock. For reliable calculations, however, a high resolution of 
the instability region is indispensable. Otherwise, the physi- 
cal strange mode instability generating acoustic energy and 
shock waves is suppressed. Thus, due to insufficient resolu- 
tion, the simulation discussed had to be stopped after ~100 
days of simulated time. To overcome the difficulty grid re- 
construction is necessary. A corresponding algorithm, con- 
sistently inserting and eliminating gridpoints during the cal- 
culation is currently being developed and will be commented 
on in a forthcoming publication. 



5 CONCLUSIONS 

Motivated by the discovery of high frequency shock oscilla- 
tions in pulsating stars confined in a narrow region, we have 
developed a procedure for an efficient treatment of such phe- 
nomena. The integration domain is decomposed into sev- 
eral subdomains according to the various, vastly different 
timescales present in the configuration. To save computing 
time, these domains are integrated according to their respec- 
tive timescales. 

Criteria for the choice and decomposition of the compu- 
tational domains have been derived. Decomposition implies 
artificial inner boundaries which require suitable boundary 
conditions. How these have to be chosen in order to minimise 
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the numerical error (compared to the approach without de- 
composition) has been discussed. An overlap of the domains 
was found to be necessary to produce satisfactory results. 

The decomposition technique has been tested and vali- 
dated by comparison with results obtained without decom- 
position. The major effect of decomposition was found to 
consist of a delay in time (phaseshift) with respect to the 
original calculation. The latter is regarded as physically 
largely irrelevant. Otherwise the numerical quality of the 
results was proven to be satisfactory. 

For the models considered, decomposition was found 
both theoretically and by numerical tests to reduce the 
computational costs by approximately a factor of two. The 
speed-up critically depends on the size of the rapidly varying 
domain. Although a reduction of computing time by 50 per 
cent sounds moderate, it is of practical relevance consider- 
ing the duration of several weeks for a complete simulation. 
The intended extension of the procedure to two dimensional 
problems will yield an appreciably higher speed-up than that 
for the one dimensional model considered here, since the iter- 
ative inversion of some matrices there requires of the order 
of iV^ operations. Moreover, decomposition in two dimen- 
sions may reduce the size of matrices to be inverted such 
that fast direct methods (requiring of the order of N opera- 
tions) become most efficient, whose application would then 
imply a further acceleration of the computation. 

We have applied the method to a model for the star P 
Cygni, paying special attention to the adequate treatment 
of the different timescales involved. 
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